Modeling the Joint Distribution of Wind Speed and Direction using Gaussain Mixture Modles¶
OEN Method: Harris, Cook The parent wind speed distribution: Why Weibull? http://www.sciencedirect.com/science/article/pii/S0167610514001056
Gaussian Mixture Models, http://scikit-learn.org/stable/modules/mixture.html
1. Set up¶
1.1 Environment¶
In [1]:
%matplotlib inline
%load_ext autoreload
%autoreload 2
# from __future__ import division
from import_file import *
from helpers.parallel_helper import *
load_libs()
1.2 Read Data¶
In [2]:
# ipath= "./data/NCDC/us/baltimore/dat.txt" # time too short
# ipath= './data/NCDC/uk/boscombe_down/dat.txt' # 4?
# ipath= './data/NCDC/uk/middle_wallop/dat.txt' # good
# ipath= './data/NCDC/uk/southhamption/dat.txt' # high 0
# ipath= './data/NCDC/uk/bournemouth/dat.txt' 4?
ipath= './data/NCDC/uk/tiree/dat.txt' # try 4
# ipath= "./data/NCDC/uk/bealach_na_ba/dat.txt" #
# ipath= "./data/NCDC/uk/skye_lusa/dat.txt" #
# ipath= "./data/NCDC/uk/south_uist_range/dat.txt"
# ipath= './data/NCDC/uk/marham/dat.txt'
# ipath= "./data/NCDC/uk/holbeach/dat.txt"
# ipath= "./data/NCDC/uk/weybourne/dat.txt"
# ipath= "./data/NCDC/uk/cambridge/dat.txt"
# ipath= "./data/NCDC/uk/wattisham/dat.txt"
# ipath= "./data/NCDC/uk/benbecula/dat.txt" # truncate (untruncate in m/s), 4?
# ipath= "./data/NCDC/southeast_asia/singapore_changi/dat.txt" # 4 or 5?
# ipath= "./data/NCDC/southeast_asia/sultan_ismail/dat.txt" #
# ipath= "./data/NCDC/southeast_asia/paya_lebar/dat.txt" # 2010-2014 has problem
# ipath= "./data/NCDC/east_asia/cheongju_intl/dat.txt" # 2005-2009 may have problem, fit is good
# ipath= "./data/NCDC/east_asia/daegu_ab/dat.txt" # recent 5 year may have problem, but fit is generally good
# ipath= "./data/NCDC/canada/melita_man/dat.txt" #
# ipath= "./data/NCDC/canada/mildred_lake_alta/dat.txt" # uptail wind speed
# ipath= "./data/NCDC/canada/manouane_est_que/dat.txt"
# ipath= "./data/NCDC/canada/argentia/dat.txt"
# ipath= "./data/NCDC/europe/landsberg_lech/dat.txt" # very good, can try 4
# ipath= "./data/NCDC/europe/vatry/dat.txt" # double peak, initial speed (should be good with m/s)
# ipath= "./data/NCDC/europe/avord/dat.txt" # try 4, initial speed (should be good with m/s)
# ipath= './data/NCDC/europe/ciampino/dat.txt' # try 4, bandwidth?
# ipath= "./data/NCDC/europe/holzdorf/dat.txt"
# ipath= "./data/NCDC/europe/neuburg/dat.txt"
# ipath= "./data/NCDC/europe/valladolid/dat.txt"
# ipath= "./data/NCDC/europe/laupheim/dat.txt" # double peak, 4; very good
# ipath= "./data/NCDC/europe/barayas/dat.txt" #
# ipath= "./data/NCDC/europe/huspel_aws/dat.txt" # integer, 4?
# ipath= "./data/NCDC/europe/diepholz/dat.txt" # recent year is integer with less data, use m/s
# ipath= './data/NCDC/europe/tenerife_sur/dat.txt' # some directions are blocked
# ipath= "./data/NCDC/oceania/narrandera/dat.txt" # high 0
# ipath= "./data/NCDC/oceania/brisbane_archerfield/dat.txt" # high 0
# ipath= "./data/NCDC/oceania/auckland_intl/dat.txt" # Weird KDE shape, 4?
# ipath= "./data/NCDC/oceania/canberra/dat.txt" # high 0
# ipath= "./data/NCDC/cn/shanghai/hongqiao_intl/dat.txt"
# ipath= "./data/NCDC/cn/shanghai/pudong/dat.txt"
# ipath= "./data/NCDC/cn/tianjing/binhai/dat.txt" # recent 5 year has small problem, good overall, few 0
# ipath= "./data/NCDC/cn/tianjing/tianjing/dat.txt" # 16 sectors
# ipath= "./data/NCDC/cn/hefei_luogang/dat.txt" # recent 5 year has small problem, good overall, few 0
# ipath= "./data/NCDC/cn/南京lukou/dat.txt"
# ipath= "./data/NCDC/cn/shijiazhuang_zhengding/dat.txt"
# ipath= "./data/NCDC/cn/henan_gushi/dat.txt" # 16 sectors, fit not very good
# ipath= "./data/NCDC/cn/郑州xinzheng/dat.txt" # fit is not good
# ipath= "./data/NCDC/cn/nanning_wuxu/dat.txt" # recent 5 year has problem
# ipath= './data/NCDC/cn/macau/dat.txt'
# ipath= "./data/NCDC/cn/hk_intl/dat.txt" # few initial data
# ipath= "./data/NCDC/us/sheboygan/dat.txt" # untruncate us dataset, but recent 5 year maybe corrupted
# ipath= "./data/cn_database/广州.txt"
# ipath= "./data/cn_database/北京.txt"
# uipath = unicode(ipath , "utf8")
uipath = ipath
if not "cn_database" in ipath:
df = pd.read_csv(uipath,header=0, skipinitialspace=True)
df.rename(columns={'Date':'date','Dir':'dir','Spd':'speed','Type':'type','I.1':'wind_type'}, inplace=True)
df = df[['date','HrMn','type','dir','speed','wind_type' ]]
else:
df = read_cn_database(uipath)
In [3]:
df = df.query("(dir <= 999) & (speed < 100) & \
(date >= 19700000) & (date < 20160000) ")
In [4]:
plot_speed_and_angle_distribution(df.speed, df.dir)
In [5]:
# Dir [10,360]=> [0,350]
df['dir'] = df['dir'].apply(lambda x: x%360 if x < 999 else x)
df['month'] = df['date']%10000//100
# Convert Windrose coordianates to Polar Cooridinates
df['dir_windrose'] = df['dir']
df['dir'] = df['dir'].apply(lambda x: (90 - x)%360 if x < 999 else x)
df.describe()
Out[5]:
In [6]:
df.plot(x='date', y='speed',legend=True,figsize=(10,3))
Out[6]:
1.2.1 Unit Detection¶
In [7]:
df['decimal'] = df.speed % 1
df.decimal.hist(alpha=0.5, label='m/s', figsize=(4, 3))
knot_unit = True if len(df.query('decimal >= 0.2')) / len(df) > 0.3 else False
# knot_unit = False
if knot_unit:
df['speed'] = df['speed'] * 1.943845
df['decimal'] = df.speed % 1
df.decimal.hist(alpha=0.5, label='knot')
# need more elaboration, some is not near an integer
df['speed'] = df['speed'].apply(lambda x: int(round(x)))
plt_configure(xlabel='Decimal', ylabel='Frequency', legend={'loc': 'best'}, title='Decimal Distribution')
print(knot_unit)
1.2.2 Sampling Type Selection¶
In [8]:
sample_type = df.query('date > 20000000')['type']
sample_type.value_counts().plot(
kind = 'bar', title = 'Report Types Comprisement', figsize=(4,3))
report_type_most_used = sample_type.value_counts().argmax()
df = df.query("type==@report_type_most_used")
1.2.3 Sampling Time Selection¶
In [9]:
df['sample_time'] = df.HrMn % 100
sample_time = df.query('date > 20000000')['sample_time']
sample_time.value_counts().plot(
kind = 'bar', title = 'sampling time distribution, in hour', figsize=(4,3))
# Need more elaboration
sample_times = sample_time.value_counts()[sample_time.value_counts() > 2000]
sample_times = sample_times.index.tolist()
# sample_times = [sample_times.index[0]]
df = df.query("sample_time in @sample_times")
df.drop(['sample_time'], 1,inplace=True)
In [10]:
df['HrMn'].value_counts().sort_index()
Out[10]:
In [11]:
MID_YEAR = (min(df.date)//10000+max(df.date)//10000)//2
bins=arange(0, 2400+100, 50)
df['HrMn'].plot(kind='hist', alpha=0.5,bins=bins, label='Overall')
df.query('date > @MID_YEAR * 10000')['HrMn'].plot(
kind='hist', alpha=0.5,bins=bins, label='> %s' % MID_YEAR )
plt_configure(xlabel='Sampling Time', ylabel='Frequency', legend={'loc':'best'}, figsize=(8, 3),
title = 'Sampling Time Distribution, Overall and > %s ' % MID_YEAR)
1.3 Data Wrangling¶
1.3.1 Artefacts¶
1.3.1.1 wrong direction record¶
In [12]:
display(df.query("(dir % 10 >= 0.1) & (dir != 999)"))
df = df.query('(dir % 10 <= 0.1) | (dir == 999)')
1.3.1.2 sudden increase in speed¶
In [13]:
# sudden increse
df['incre'] = df.speed.diff(1)
df['incre'].fillna(0, inplace=True)
df['incre_reverse'] = df.speed.diff(-1)
df['incre_reverse'].fillna(0, inplace=True)
display(df.sort_values(by='speed',ascending=False).head(10))
df['incre'].plot(kind='hist', bins=arange(-15, 15), legend=True, figsize=(8, 3))
Out[13]:
In [14]:
incre_threshold = 20 if knot_unit else 10
print('sudden increase number', len(df.query('(incre > @incre_threshold )&(incre_reverse > @incre_threshold )')))
df = df.query('(incre < @incre_threshold )|(incre_reverse < @incre_threshold )')
After Simple Aretefacts Handling
In [15]:
# Check the max speed
df.sort_values(by='speed',ascending=False).head(10)
Out[15]:
1.3.2 Direction re-aligment and 999¶
For some dataset, the 16 sectors are not record properly,
e.g. the sectors are [0,20,30,50], need to redistribute the angle into 22.5
In [16]:
display(df['dir'].value_counts().sort_index())
effective_column = df.query('dir < 999')['dir'].value_counts()[df['dir'].value_counts() > 30].sort_index()
SECTOR_LENGTH = 360/len(effective_column)
print(len(effective_column), SECTOR_LENGTH)
In [17]:
df=realign_direction(df, effective_column)
In [18]:
df=fill_direction_999(df, SECTOR_LENGTH)
1.3.3 Time Shift Comparison¶
In [19]:
DIR_REDISTRIBUTE = 'even'
if DIR_REDISTRIBUTE == 'even':
DIR_BIN = arange(-5, 360, 10)
elif DIR_REDISTRIBUTE == 'round_up':
DIR_BIN = arange(0, 360+10, 10)
# Comparison between mid_year, looking for:
# 1. Odd Even Bias
# 2. Time Shift of Wind Speed Distribution
bins = arange(0, max(df.speed) + 1)
df.query('date < @MID_YEAR * 10000')['speed'].plot(
kind='hist', alpha=0.5,bins=bins, label='< %s' % MID_YEAR)
df.query('date > @MID_YEAR * 10000')['speed'].plot(
kind='hist', alpha=0.5,bins=bins, label='> %s' % MID_YEAR)
plt.suptitle('Speed Comparison between year < %s, > %s ' % (MID_YEAR, MID_YEAR), fontsize = 14)
plt_configure(xlabel='Speed', ylabel='Frequency', legend=True, figsize=(8, 3))
In [20]:
df.query('date < @MID_YEAR * 10000')['dir'].plot(
kind='hist', alpha=0.5,bins=DIR_BIN, label='< %s' % MID_YEAR)
df.query('date > @MID_YEAR * 10000')['dir'].plot(
kind='hist', alpha=0.5,bins=DIR_BIN, label='> %s' % MID_YEAR)
plt.suptitle('Dir Comparison between year < %s, and > %s ' % (MID_YEAR, MID_YEAR), fontsize = 14)
plt_configure(xlabel='Dir', ylabel='Frequency', legend={'loc':'best'}, figsize=(8, 3),tight='x')
In [21]:
# Inspect the time shift of speed and degree distribution, and odd-even bias
check_time_shift(df)
1.3.4 0 Speed¶
In [22]:
with_too_many_zero, null_wind_frequency = is_with_too_many_zero(df.query("(date >= 20050000)"))
delete_zero = with_too_many_zero
if delete_zero:
df = df.query('(speed > 0)')
print(delete_zero, null_wind_frequency)
1.4 Re-distribute Direction and Speed¶
e.g. Dir 50 -> -45 ~ 55, to make KDE result better
In [23]:
# Randomize the angle degree, to form a better KDE estimation
df = randomize_angle(df, DIR_REDISTRIBUTE, SECTOR_LENGTH)
In [24]:
if delete_zero:
redistribute_method = 'round_down'
else:
redistribute_method = 'round_up'
df, speed_redistribution_info = randomize_speed(df, redistribute_method)
1.5 Generate (x,y) from (speed,dir)¶
In [25]:
# Cook orientation
# df['dir']= (df['dir'] + 180)%360
In [26]:
# There might be a small dot in the centre, which is due to too many zero (more than 1 speed) in center
# Scatter plot in matplot has performance issue, the speed is very slow
df['x'] = df['speed'] * cos(df['dir'] * pi / 180.0)
df['y'] = df['speed'] * sin(df['dir'] * pi / 180.0)
2. Re-select Data and Configs¶
2.1 Data Overview¶
In [27]:
## Summery of the data selection
print('Knot unit?', knot_unit)
print('Incre_threshold for artefact detection:', incre_threshold)
print('Report type used:', report_type_most_used)
print('Sampling time used:', sample_times)
print('delete_zero?', delete_zero)
if 'speed_redistribution_info' in globals():
print('Speed redistribution info:', speed_redistribution_info )
df_all_years = df # for later across-year comparison
df = df_all_years.query('(date >= 20100000) & (date < 20150000)')
# df = df.query('(HrMn == 0) and (speed >= 0.5) and (date%10000 > 900) and (date%10000 < 1000)' )
df.describe()
Out[27]:
In [28]:
# 90 degree is east
ax = WindroseAxes.from_ax()
viridis = plt.get_cmap('viridis')
ax.bar(df.dir_windrose, df.speed, normed=True, opening=0.8, edgecolor='white', nsector=36, cmap=viridis)
ax.set_legend()
In [29]:
ax, _params = WindAxes.from_ax().pdf(df.speed, bins=arange(0, max(df.speed)))
weibull_params = sp.stats.exponweib.fit(df.speed, floc=0, f0=1)
plt_configure(figsize=(4,3),xlabel='V',ylabel='PDF')
fig,ax = plt.subplots()
ecdf = sm.distributions.ECDF(df.speed)
x = linspace(0, df.speed.max())
y = ecdf(x)
plot(log(x), log(-log(1-y)),'o', label='ECDF')
k_shape, lamb_scale = weibull_params[1], weibull_params[3]
y=1-exp(-(x/lamb_scale)**k_shape) # cdf
plot(log(x), log(-log(1-y)),'-', label='Weibull')
plt_configure(xlabel="ln(V)", ylabel="ln(-ln(1-P)",legend={'loc':'best'},figsize=(4,3))
align_figures()
In [30]:
df.plot(kind='scatter', x='x', y='y', alpha=0.10, s=2)
plt.gca().set_aspect('equal')
2.2 Overview by Month¶
In [31]:
viridis = plt.get_cmap('viridis')
month_incre = 1
current_df = df.query('speed>=1')
for month in arange(1, 12+month_incre, month_incre):
end_month = month+month_incre
sub_df = current_df.query('(month >= @month) and (month < @end_month)')
if len(sub_df) > 0:
if month_incre == 1:
title = 'Month: %s' % (month)
else:
title = 'Month: %s - %s ' % (month, end_month-1)
ax = WindroseAxes.from_ax()
ax.bar(sub_df.dir_windrose, sub_df.speed, normed=True, opening=0.8, edgecolor='white', nsector=36, cmap=viridis)
plt_configure(figsize=(3,3), title=title)
# fig,ax = plt.subplots(figsize=(3,3))
# sub_df['speed'].hist(bins=arange(0,sub_df.speed.max()))
align_figures()
2.3 Create input data, i.e. speed_set¶
In [32]:
SPEED_SET = array(list(zip(df.x, df.y)))
NUMBER_OF_GAUSSIAN = 3
FIT_METHOD = 'square_error'
KDE_GRIDSEARCH = True
DEFAULT_BANDWDITH = 1.5 if knot_unit else 0.7
fig_list = []
In [33]:
FITTING_RANGE = []
fit_limit = ceil(df['speed'].quantile(.95))
fitting_axis_range = arange(-fit_limit, fit_limit+1, 1)
print(fitting_axis_range)
for i in fitting_axis_range:
for j in fitting_axis_range:
FITTING_RANGE.append([i,j])
In [34]:
plot_limit = ceil(df['speed'].quantile(.95))
PLOT_AXIS_RANGE = arange(-plot_limit, plot_limit+1, 1)
3. Kernel Density Estimation¶
In [35]:
sample = SPEED_SET
KDE_KERNEL = 'gaussian'
# KDE_KERNEL, bandwidth = 'tophat', 1
In [36]:
# %%time
# from sklearn.grid_search import GridSearchCV
# if KDE_GRIDSEARCH == True:
# # The bandwidth value sometimes would be too radical
# if knot_unit:
# bandwidth_range = arange(0.7,2,0.2)
# else:
# bandwidth_range = arange(0.4,1,0.1)
# # Grid search is unable to deal with too many data (a long time is needed)
# if len(sample) > 50000:
# df_resample=df.sample(n=50000, replace=True)
# bandwidth_search_sample = array(list(zip(df_resample.x, df_resample.y)))
# else:
# bandwidth_search_sample = sample
# grid = GridSearchCV(neighbors.KernelDensity(kernel = KDE_KERNEL),
# {'bandwidth': bandwidth_range},
# n_jobs=-1, cv=4)
# grid.fit(bandwidth_search_sample)
# bandwidth = grid.best_params_['bandwidth']
In [37]:
if 'bandwidth' not in globals():
bandwidth = DEFAULT_BANDWDITH
bandwidth = 1.9
kde = neighbors.KernelDensity(bandwidth=bandwidth,kernel = KDE_KERNEL).fit(sample)
points = FITTING_RANGE
# very slow if the dataset is too large, e.g. 100,000
# kde returns log prob, need to convert it
kde_result = exp(kde.score_samples(points))
print('bandwidth:', bandwidth, len(kde_result))
print(kde_result[:5])
In [38]:
# Plot jPDF
X = Y = PLOT_AXIS_RANGE
# Can't work if pass as generate_Z_from_X_Y(X,Y, exp(kde.score_samples())), need to use lambda
# see http://stackoverflow.com/questions/21035437/passing-a-function-as-an-argument-in-python
kde_Z = generate_Z_from_X_Y(X,Y, lambda coords: exp(kde.score_samples(coords)))
colorbar_lim = 0, kde_Z.max()
plot_3d_prob_density(X,Y,kde_Z)
fig_kde,ax1 = plt.subplots(figsize=(4,3))
plot_2d_prob_density(X,Y,kde_Z,xlabel='x', ylabel='y',ax=ax1)
with sns.axes_style({'axes.grid' : False}):
from matplotlib import ticker
fig_hist,ax2 = plt.subplots(figsize=(4,3))
_,_,_,image = ax2.hist2d(df.x, df.y, bins=PLOT_AXIS_RANGE, cmap='viridis',)
ax2.set_aspect('equal')
cb = plt.colorbar(image)
tick_locator = ticker.MaxNLocator(nbins=6)
cb.locator = tick_locator
cb.update_ticks()
plt_configure(ax=ax2, xlabel='x',ylabel='y')
align_figures()
In [39]:
kde_cdf = cdf_from_pdf(kde_result)
4. GMM by Expectation-maximization¶
In [40]:
sample= SPEED_SET
clf = mixture.GMM(n_components=NUMBER_OF_GAUSSIAN, covariance_type='full')
clf.fit(sample)
print(clf.converged_)
In [41]:
gmm_em_result = read_gmm_em_result(clf)
pretty_print_gmm(gmm_em_result)
Out[41]:
In [42]:
fig,ax = plt.subplots(figsize=(3.5,3.5))
plot_gmm_ellipses(gmm_em_result, ax=ax)
In [43]:
X = Y = PLOT_AXIS_RANGE
pdf_Z = generate_Z_from_X_Y(X,Y, lambda coords: exp(clf.score_samples(coords)[0]))
def residule_between_kde_and_gmm(points):
kde_vals = exp(kde.score_samples(points))
gmm_vals = exp(clf.score_samples(points)[0])
return kde_vals - gmm_vals
residual_Z = generate_Z_from_X_Y(X,Y, residule_between_kde_and_gmm)
plot_3d_prob_density(X,Y,pdf_Z)
plot_3d_prob_density(X,Y,residual_Z)
align_figures()
fig = plt.figure(figsize=(4,3))
plot_2d_prob_density(X,Y,kde_Z,xlabel='x',ylabel='y', colorbar_lim=colorbar_lim)
fig_em = plt.figure(figsize=(4,3))
plot_2d_prob_density(X,Y,pdf_Z,xlabel='x',ylabel='y', colorbar_lim=colorbar_lim)
fig = plt.figure(figsize=(4,3))
plot_2d_prob_density(X,Y,residual_Z,xlabel='x',ylabel='y')
align_figures()
Goodness-of-fit Statistics¶
In [44]:
points = FITTING_RANGE
gmm_pdf_result = exp(clf.score_samples(points)[0])
gof_df(gmm_pdf_result, kde_result)
Out[44]:
5. GMM by Optimization¶
In [45]:
sample = SPEED_SET
points = FITTING_RANGE
max_speed = df.speed.max()
print(FIT_METHOD)
In [46]:
# from GMM,EM
# GMM format: weight, meanx, meany, sigx, sigy, rho
x0 = gmm_em_result
cons = [
# sum of every 6th element, which is the fraction of each gaussian
{'type': 'eq', 'fun': lambda x: sum(x[::6]) - 1},
# # limit the width/height ratio of elliplse, optional
# {'type': 'ineq', 'fun': lambda x: width_height_ratios_set(x) - 1/3},
# {'type': 'ineq', 'fun': lambda x: 3 - width_height_ratios_set(x)},
]
bonds = [(0., 0.99),(-fit_limit, fit_limit),(-fit_limit, fit_limit),
(0., fit_limit),(0., fit_limit),(-0.99, 0.99)]*(len(x0)//6)
result = sp.optimize.minimize(
lambda x0: GMM_fit_score(x0, kde_result, points, FIT_METHOD),
x0,
bounds = bonds,
constraints=cons,
tol = 0.000000000001,
options = {"maxiter": 500})
result
Out[46]:
5.1 GMM Result¶
In [47]:
gmm = group_gmm_param_from_gmm_param_array(result.x, sort_group = True)
mixed_model_pdf = generate_gmm_pdf_from_grouped_gmm_param(gmm)
gmm_pdf_result = mixed_model_pdf(points)
pretty_print_gmm(gmm)
Out[47]:
In [48]:
fig_gmm, ax = plt.subplots(figsize=(3.5,3.5))
plot_gmm_ellipses(gmm,ax=ax)
# fig_gmm.savefig('foo.pdf')
5.2 Goodness-of-fit statistics¶
In [49]:
gof_df(gmm_pdf_result, kde_result)
Out[49]:
In [50]:
X = Y = PLOT_AXIS_RANGE
pdf_Z = generate_Z_from_X_Y(X,Y, mixed_model_pdf)# passing a function as an argument
def residule_between_kde_and_gmm(points):
kde_vals = exp(kde.score_samples(points))
gmm_vals = mixed_model_pdf(points)
return kde_vals - gmm_vals
residual_Z = generate_Z_from_X_Y(X,Y, residule_between_kde_and_gmm)
plot_3d_prob_density(X,Y,pdf_Z)
plot_3d_prob_density(X,Y,residual_Z)
align_figures()
fig = plt.figure(figsize=(4,3))
plot_2d_prob_density(X,Y,kde_Z,xlabel='x',ylabel='y')
fig_gmm = plt.figure(figsize=(4,3))
plot_2d_prob_density(X,Y,pdf_Z,xlabel='x',ylabel='y')
fig = plt.figure(figsize=(4,3))
plot_2d_prob_density(X,Y,residual_Z,xlabel='x',ylabel='y')
align_figures()
In [51]:
def f(V,theta):
return (mixed_model_pdf([[V*cos(theta),V*sin(theta)]]))*V
In [52]:
x = linspace(0, max_speed, 15)
# 1. Plot v-p points, empirical distribution
ecdf = sm.distributions.ECDF(df.speed)
y = ecdf(x)
plot(log(x), log(-log(1-y)),'o', label = 'Empirical')
# 2. Weibull distribution
weibull_params = sp.stats.exponweib.fit(df.speed, floc=0, f0=1)
y_cdf_weibull = sp.stats.exponweib.cdf(x, *weibull_params)
plot(log(x), log(-log(1-y_cdf_weibull)),'--', label = 'Weibull')
# 3. GMM distribution
y_ = [integrate.nquad(f, [[0, x_val],[0, 2*pi]]) for x_val in x]
y = array(list(zip(*y_))[0])
plot(log(x), log(-log(1-y)),'-', color='black', label = 'GMM')
plt_configure(xlabel='ln(V)',ylabel='ln(-ln(1-P))',legend={'loc':'best'})
In [53]:
# Calculate Angle Distribution
x = linspace(0,2*pi, num=36+1)
y_ =[integrate.nquad(f, [[0, inf],[x_-pi/36, x_+pi/36]]) for x_ in x]
# 1. integrate.nquad returns 2 vaule, result, abserr. Need to abstract the first to plot
# 2. * length of data size, to adjust the plot space
y = array(list(zip(*y_))[0])*len(df['dir'])
df['dir'].hist(bins=DIR_BIN, alpha=0.5, label='Data')
plot(x/pi*180, y,'-', color='black', label='GMM')
title='Direction Distribution Comparison'
plt_configure(xlabel='Direction',ylabel='Frequency',
legend={'loc': 'best'} ,tight='xtight',figsize = (4.5,2))
dir_fig = plt.gcf()
print(title)
In [54]:
# Calculate Speed Distribution
# 1. GMM Model
x = arange(0, max_speed, 0.5)
y_ =[integrate.nquad(f, [[x_-0.01, x_+0.01],[0, 2*pi]]) for x_ in x]
y = array(list(zip(*y_))[0])*len(df.speed)/0.02
plot(x, y,'-', color='black', label='GMM')
# 2. Observation
df['speed'].hist(bins=arange(0, df.speed.max()), alpha=0.5, label='Data')
# 3. Weibull
y_weibul = sp.stats.exponweib.pdf(x, *weibull_params)* len(df.speed)
plot(x, y_weibul, '--', color='black', label='Weibull')
title='Speed Distribution Comparison'
plt_configure(xlabel='Speed',ylabel='Frequency',legend=True, figsize=(4.5, 2))
plt.gca().set_ylim(bottom = 0)
plt.locator_params(axis='y', nbins=5)
print(title)
In [55]:
# %%time
incre = SECTOR_LENGTH
density_collection=Parallel(n_jobs=-1)(delayed(direction_compare)(gmm, df, angle, incre)
for angle in arange(0, 360, incre))
# #This R square is compute as in paper
# # Comparison of bivariate distribution constructionapproaches for analysing wind speed anddirection data
# # http://onlinelibrary.wiley.com/doi/10.1002/we.400/full
print(true_R_square(density_collection))
5.3 Sectoral Comaprison¶
In [56]:
# Calculate Speed Distribution
def model_data_comparison(df, original_incre = 10, incre = 10):
start, end = -original_incre/2 + incre/2, 360
max_diff_array = []
curve_collection = []
max_speed = df.speed.max()
# Find a max count for plotting histogram
max_count = max_count_for_angles(df, start, end, incre)
plot_range = [0, max_speed, 0, max_count*1.05]
for angle in arange(start, end, incre):
angle_radian, incre_radian = radians(angle), radians(incre)
start_angle, end_angle = angle-incre/2, angle+incre/2
# Select data from observation
sub_df, sub_max_speed = select_df_by_angle(df, start_angle, end_angle)
data_size = len(sub_df.speed)
direction_prob = integrate.nquad(f, [[0, inf],
[angle_radian-incre_radian/2, angle_radian+incre_radian/2]])
weibull_params = sp.stats.exponweib.fit(sub_df.speed, floc=0, f0=1)
fig = plt.figure(figsize=(13,1.3))
x = linspace(0, sub_max_speed, 20)
# 1. Frequency Comparison
ax1 = fig.add_subplot(1,3,1)
# 1.1. Data Histogram
sub_df['speed'].hist(bins=arange(0, sub_max_speed), alpha=0.5, label='Data')
# 1.2. Model Distribution
y_gmm_ =[integrate.nquad(f, [[x_-0.01, x_+0.01],
[angle_radian-incre_radian/2, angle_radian+incre_radian/2]])
for x_ in x]
y_gmm = array(list(zip(*y_gmm_))[0])/direction_prob[0]/0.02
plot(x, y_gmm*data_size,'-', color='black', label='GMM')
# 1.3 Weibull Distribution
y_weibull = sp.stats.exponweib.pdf(x, *weibull_params)
plot(x, y_weibull*data_size, '--', color='black',label='Weibull')
plt_configure(xlabel = "V", ylabel='Frequency', legend=True)
plt.axis(plot_range)
# 2. CDF Comaprison
ax2 = fig.add_subplot(1,3,2)
# 2.1. Empirical distribution
ecdf = sm.distributions.ECDF(sub_df.speed)
y_ecdf = ecdf(x)
# 2.2. Model
y_ =[integrate.nquad(f, [[0, x_val],[angle_radian-incre_radian/2, angle_radian+incre_radian/2]])
for x_val in x]
y_cdf_gmm = array(list(zip(*y_))[0])/direction_prob[0]
# 2.3 Weibull
y_cdf_weibull = sp.stats.exponweib.cdf(x, *weibull_params)
plot(x, y_ecdf,'o', alpha=0.8, label='Data')
plot(x, y_cdf_gmm,'-', color='black',label='GMM')
plot(x, y_cdf_weibull,'--', color='black',label='Weibull')
plt.gca().set_xlim(right = max_speed)
plt_configure(xlabel = "V", ylabel='P', legend=True)
# 3. Weibull Comparison
ax3 = fig.add_subplot(1,3,3)
plot(log(x), log(-log(1-y_ecdf)),'o', alpha=0.8, label='Data')
plot(log(x), log(-log(1-y_cdf_gmm)),'-', color='black', label='GMM')
plot(log(x), log(-log(1-y_cdf_weibull)),'--',color='black',label='Weibull')
plt_configure(xlabel="ln(V)", ylabel="ln(-ln(1-P)",legend={'loc':'best'})
plt.gca().set_xlim(right = log(max_speed+1))
print('%s (%s - %s) degree' % (angle, start_angle, end_angle) )
bins = arange(0, sub_df.speed.max()+1)
density, _ = np.histogram(sub_df['speed'],bins=bins, normed=True)
density_expected_ =[integrate.nquad(f, [[x_, x_+1],[angle_radian-incre_radian/2, angle_radian+incre_radian/2]])
for x_ in bins[:-1]]
density_expected_gmm = array(list(zip(*density_expected_ ))[0])/direction_prob[0]
R_square_gmm = sector_r_square(density, density_expected_gmm)
density_expected_weibull = sp.stats.exponweib.cdf(bins[1:], *weibull_params) - sp.stats.exponweib.cdf(bins[:-1], *weibull_params)
R_square_weibull = sector_r_square(density, density_expected_weibull)
diff, diff_weibull= np.abs(y_ecdf - y_cdf_gmm), np.abs(y_ecdf - y_cdf_weibull)
max_diff_array.append([len(sub_df), angle, diff.max(), x[diff.argmax()],
diff_weibull.max(), x[diff_weibull.argmax()], R_square_gmm, R_square_weibull])
curves = {'angle': angle, 'data_size': data_size, 'weight': direction_prob[0],
'x': x, 'gmm_pdf': y_gmm, 'gmm_cdf': y_cdf_gmm,
'weibull_pdf': y_weibull, 'weibull_cdf': y_cdf_weibull, 'ecdf': y_ecdf}
curve_collection.append(curves)
plt.show()
print('data size:', len(sub_df), 'weight', len(sub_df)/len(df))
print('GMM', 'Weibull')
print('R square', R_square_gmm, R_square_weibull)
print('max diff:', diff.max(), diff_weibull.max(), 'speed value:', x[diff.argmax()], 'y gmm', y_cdf_gmm[diff.argmax()])
print(' ')
return max_diff_array, curve_collection
In [57]:
%%time
if len(effective_column) == 16:
rebinned_angle = 22.5
else:
rebinned_angle = 20
max_diff_array, curve_collection = model_data_comparison(df, SECTOR_LENGTH, rebinned_angle)
In [58]:
diff_df = pd.DataFrame(max_diff_array,columns=['datasize','direction', 'gmm', 'speed_gmm',
'weibull', 'speed_weibull', 'r_square_gmm', 'r_square_weibull'])
gmm_mean, weibull_mean = plot_sectoral_comparison(diff_df.r_square_gmm, diff_df.r_square_weibull, diff_df.direction, diff_df.datasize)
plt_configure(ylabel="$\ R^2$")
print(gmm_mean, weibull_mean)
In [59]:
gmm_mean, weibull_mean = plot_sectoral_comparison(diff_df.gmm, diff_df.weibull, diff_df.direction, diff_df.datasize)
plt_configure(ylabel="K-S")
ylim = max(plt.gca().get_ylim()[1],0.15)
plt.gca().set_ylim(top=ylim)
print(gmm_mean, weibull_mean)
In [60]:
# Compare direction weight with previous figure
display(dir_fig)
5.3.1 Data Variability, by Bootstrap (Resampling)¶
In [61]:
max_diff_element = max(max_diff_array, key=lambda x: x[2])
angle = max_diff_angle = max_diff_element[1]
incre = rebinned_angle
In [62]:
FRACTION = 1
# Select data from observation
start_angle, end_angle = angle-incre/2, angle+incre/2
angle_radian, incre_radian = radians(angle), radians(incre)
sub_df, sub_max_speed = select_df_by_angle(df, start_angle, end_angle)
weibull_params = sp.stats.exponweib.fit(sub_df.speed, floc=0, f0=1)
In [63]:
fig = plt.figure(figsize=(13,1.5))
ax1 = fig.add_subplot(1,3,1)
ax2 = fig.add_subplot(1,3,2)
ax3 = fig.add_subplot(1,3,3)
x = arange(0, sub_max_speed, 0.5)
# 1. Data
bins=arange(0, sub_max_speed)
sub_df['speed'].hist(ax=ax1, bins=bins, alpha=0.5, label='Data', normed=True)
# 2. GMM
y_gmm_ =[integrate.nquad(f, [[x_-0.01, x_+0.01], [angle_radian-incre_radian/2, angle_radian+incre_radian/2]])
for x_ in x]
direction_prob = integrate.nquad(f, [[0, inf], [angle_radian-incre_radian/2, angle_radian+incre_radian/2]])
y_gmm = array(list(zip(*y_gmm_))[0])/direction_prob[0]/0.02
y_gmm_cdf_ =[integrate.nquad(f, [[0, x_val],[angle_radian-incre_radian/2, angle_radian+incre_radian/2]]) for x_val in x]
y_gmm_cdf = array(list(zip(*y_gmm_cdf_))[0])/direction_prob[0]
ax1.plot(x, y_gmm,'-', color='black', label='GMM')
ax2.plot(x, y_gmm_cdf,'-', color = 'black', label='GMM')
ax3.plot(log(x), log(-log(1-y_gmm_cdf)),'-', color = 'black',label='GMM')
# 3. Weilbull
y_weibull = sp.stats.exponweib.pdf(x, *weibull_params)
y_weibull_cdf = sp.stats.exponweib.cdf(x, *weibull_params)
ax1.plot(x, y_weibull,'--',color='black',label='Weibull')
ax2.plot(x, y_weibull_cdf,'--',label='Weibull')
ax3.plot(log(x), log(-log(1-y_weibull_cdf)),'--',label='Weibull')
# 4. Data Resampled
count_collection = []
for i in range(1,100):
sub_df_resampled = sub_df.sample(frac=FRACTION, replace=True)
resampled_count, _ = np.histogram(sub_df_resampled['speed'], bins=bins, normed=True)
count_collection.append(resampled_count)
ecdf = sm.distributions.ECDF(sub_df_resampled.speed)
y_ecdf = ecdf(x)
ax2.plot(x, y_ecdf,':', label='Data Resampled')
ax3.plot(log(x), log(-log(1-y_ecdf)),':', label='Data Resampled')
if i == 1:
plt_configure(ax=ax2, xlabel='V',ylabel='P',legend={'loc':'best'})
plt_configure(ax=ax3, xlabel='V',ylabel='ln(-ln(1-P)',legend={'loc':'best'})
print('%s (%s - %s) Degree Speed Distribution' % (angle, start_angle, end_angle))
count_collection = np.array(count_collection)
mx, mn = np.max(count_collection,0), np.min(count_collection,0)
ax1.plot(bins[1:]-0.5, mx, ':', color='blue')
ax1.plot(bins[1:]-0.5, mn, ':', color='blue', label='Resample limit')
ax1.set_ylim(bottom = 0)
plt_configure(ax=ax1, xlabel='V',ylabel='Frequency',legend={'loc':'best'})
ax1.locator_params(axis='y', nbins=5)
ax2.locator_params(axis='y', nbins=5)
ax3.locator_params(axis='y', nbins=5)
diff = abs(y_ecdf - y_gmm_cdf)
print(diff.max(), x[diff.argmax()], y_gmm_cdf[diff.argmax()])
5.3.2 Time Variability¶
In [64]:
fig_time_variability_3d = plt.figure()
ax1 = fig_time_variability_3d.gca(projection='3d')
fig_time_variability_cdf,ax2 = plt.subplots(figsize=(3,1.8))
fig_time_variability_weibull, ax3 = plt.subplots(figsize=(3,1.8))
x = arange(0, sub_max_speed, 0.5)
# 1. GMM
direction_prob = integrate.nquad(f, [[0, inf], [angle_radian-incre_radian/2, angle_radian+incre_radian/2]])
y_gmm_ =[integrate.nquad(f, [[x_-0.01, x_+0.01], [angle_radian-incre_radian/2, angle_radian+incre_radian/2]])
for x_ in x]
y_gmm = array(list(zip(*y_gmm_))[0])/direction_prob[0]/0.02
y_gmm_cdf_ =[integrate.nquad(f, [[0, x_val],
[angle_radian-incre_radian/2, angle_radian+incre_radian/2]]) for x_val in x]
y_gmm_cdf = array(list(zip(*y_gmm_cdf_))[0])/direction_prob[0]
ax2.plot(x, y_gmm_cdf,'-', color='black', label = 'GMM')
ax3.plot(log(x), log(-log(1-y_gmm_cdf)),'-', color='black',label='GMM')
# 2. Weibull
y_weibull = sp.stats.exponweib.pdf(x, *weibull_params)
y_weibull_cdf = sp.stats.exponweib.cdf(x, *weibull_params)
ax2.plot(x, y_weibull_cdf,'--', label='Weibull')
ax3.plot(log(x), log(-log(1-y_weibull_cdf)), '--', label='Weibull')
# 3. Data
prop_cycle=iter(mpl.rcParams['axes.color_cycle'])
for start_time in range(20000000, 20150000, 50000):
end_time = start_time + 50000
time_label = start_time//10000
df_other_years = df_all_years.query('(date >= @start_time) & (date < @end_time)')
df_other_years_at_angle, sub_max_speed_other_year = select_df_by_angle(df_other_years, start_angle, end_angle)
if len(df_other_years_at_angle) > 0 :
ecdf = sm.distributions.ECDF(df_other_years_at_angle.speed)
y_ecdf = ecdf(x)
ax2.plot(x, y_ecdf,':', label = time_label)
ax3.plot(log(x), log(-log(1-y_ecdf)),':', label = time_label)
title = '%s - %s' %(time_label, time_label+4)
count, division = np.histogram(df_other_years_at_angle['speed'], normed=True,
bins=arange(0, sub_max_speed_other_year))
ax1.bar(left=division[:-1], height=count, zs=time_label, zdir='x',
color=next(prop_cycle), alpha=0.8)
x_3d = time_label*np.ones_like(x)
ax1.plot(x_3d, x, y_gmm, '-', color='black', label='GMM' if time_label == 2010 else '')
ax1.plot(x_3d, x, y_weibull, '--', color='blue', label='Weibull' if time_label == 2010 else '')
print('%s (%s - %s) Degree Speed Distribution' % (angle, start_angle, end_angle))
ax1.set_ylim(bottom = 0)
ax1.set_zlabel('Frequency')
plt_configure(ax=ax1, xlabel='Time',ylabel='V', legend=True)
plt_configure(ax=ax2, xlabel='V',ylabel='P',legend={'loc':'best'})
plt_configure(ax=ax3, xlabel="ln(V)", ylabel="ln(-ln(1-P)",legend={'loc':'best'})
ax1.set_zlim(bottom = 0)
align_figures()
5.3.3 Adjacent Sector Variability¶
In [65]:
incre = rebinned_angle
angle_group = [max_diff_angle-incre, max_diff_angle, max_diff_angle+incre]
In [66]:
fig_adjecent_variability_3d = plt.figure()
ax1 = fig_adjecent_variability_3d.gca(projection='3d')
fig_adjecent_variability_cdf, ax2 = plt.subplots(figsize=(3,1.8))
fig_adjecent_variability_weibull, ax3 = plt.subplots(figsize=(3,1.8))
legend_3d = False
prop_cycle=iter(mpl.rcParams['axes.color_cycle'])
curve_df = pd.DataFrame(curve_collection)
for angle in angle_group:
curves = curve_df.query('angle == @angle%360').T.to_dict()
curves = curves[list(curves)[0]]
data_size, x = curves['data_size'], curves['x']
y_gmm, y_cdf_gmm = curves['gmm_pdf'], curves['gmm_cdf']
y_weibull, y_cdf_weibull, y_cdf = curves['weibull_pdf'], curves['weibull_cdf'], curves['ecdf']
linestyle = '-' if angle == max_diff_angle else ':'
alpha = 0.7 if angle == max_diff_angle else 0.3
ax2.plot(x, y_gmm*data_size, linestyle, label=angle)
ax3.plot(x, y_weibull*data_size, linestyle, label=angle)
start_angle, end_angle = angle-incre/2, angle+incre/2
sub_df, sub_max_speed = select_df_by_angle(df, start_angle, end_angle)
x_3d = angle*np.ones_like(x)
ax1.plot(x_3d, x, y_gmm*data_size, color='black', label='GMM')
ax1.plot(x_3d, x, y_weibull*data_size, color='blue', linestyle='--',label='Weibull')
count, division = np.histogram(sub_df['speed'], bins=arange(0, sub_max_speed))
ax1.bar(left=division[:-1], height=count, zs=angle, zdir='x', color=next(prop_cycle), alpha=0.8)
if legend_3d == False:
ax1.legend()
legend_3d = True
plt_configure(ax=ax1, xlabel='Direction', ylabel='Speed')
plt_configure(ax=ax2, xlabel='V',ylabel='Frequency',legend={'loc':'best'})
plt_configure(ax=ax3, xlabel='V',ylabel='Frequency',legend={'loc':'best'})
ax1.set_zlabel('Frequency')
ax1.set_zlim(bottom = 0)
ylim = max(ax1.get_ylim()[1],ax3.get_ylim()[1])
ax2.set_ylim(bottom = 0, top=ylim)
ax3.set_ylim(bottom = 0, top=ylim)
print(max_diff_angle)
print('GMM, Weibull, Histogram')
align_figures()
6. Result Variability & Cross-Validation¶
In [67]:
if 'bandwidth' not in globals():
bandwidth = DEFAULT_BANDWDITH
if 'FIT_METHOD' not in globals():
FIT_METHOD = 'square_error'
if 'KDE_KERNEL' not in globals():
KDE_KERNEL = 'gaussian'
config = {'bandwidth': bandwidth,
'fitting_range': FITTING_RANGE,
'fit_limit': fit_limit,
'kde_kernel': KDE_KERNEL}
6.1 Variability of the Result¶
In [68]:
%%time
results = Parallel(n_jobs=-1)(delayed(resampled_fitting)(df, FIT_METHOD, NUMBER_OF_GAUSSIAN, config) for i in range(2))
for result in results:
display(pretty_print_gmm(result['gmm']))
fig,ax = plt.subplots(figsize=(3.5,3.5))
plot_gmm_ellipses(result['gmm'],ax=ax)
plt.show()
# fig_gmm = plt.figure(figsize=(4,3))
# mixed_model_pdf = generate_gmm_pdf_from_grouped_gmm_param(result['gmm'])
# pdf_Z = generate_Z_from_X_Y(X,Y, mixed_model_pdf)
# plot_2d_prob_density(X,Y,pdf_Z,xlabel='x',ylabel='y')
# align_figures()
display(gof_df(result['gmm_pdf_result'], result['kde_result']))
display(gof_df(result['gmm_pdf_result'], kde_result))
print('')
6.2 Cross-validation, to select the number of Gaussian¶
In [69]:
%%time
from sklearn.cross_validation import train_test_split, KFold
## 5-fold cross validation
gaussian_number_range = arange(1,6)
CV_result_train_all,CV_result_test_all =[],[]
number_of_fold = 4
print('Number of train/test dataset', len(df)*(number_of_fold-1)/number_of_fold, len(df)/number_of_fold)
for number_of_gaussian in gaussian_number_range:
print( ' ')
print('Number of gaussian', number_of_gaussian)
kf = KFold(len(df), n_folds=number_of_fold, shuffle=True)
CV_result = Parallel(n_jobs=-1)(delayed(fit_per_fold)(df, train_index, test_index, FIT_METHOD, number_of_gaussian, config) for train_index, test_index in kf)
CV_result_train, CV_result_test = list(zip(*CV_result))
CV_result_train, CV_result_test = list(CV_result_train), list(CV_result_test)
CV_result_train_all.append(CV_result_train)
CV_result_test_all.append(CV_result_test)
print('Train')
pretty_pd_display(CV_result_train)
print('Test')
pretty_pd_display(CV_result_test)
In [70]:
train_scores_mean, train_scores_std = generate_mean_std_gof(CV_result_train_all)
print('Train gof mean, std')
display(train_scores_mean)
test_scores_mean, test_scores_std = generate_mean_std_gof(CV_result_test_all)
print('Test gof mean, std')
display(test_scores_mean)
In [71]:
prop_cycle=mpl.rcParams['axes.color_cycle']
gaussian_number_range = train_scores_mean.index
for column, column_name in zip(['R_square','K_S','Chi_square'],["$\ R^2$", "K-S", "$\widetilde{\chi^2} $"]):
plot(gaussian_number_range, train_scores_mean[column],
'-', label = 'train', color=prop_cycle[0])
plt.fill_between(gaussian_number_range,
train_scores_mean[column] - train_scores_std[column],
train_scores_mean[column] + train_scores_std[column],
alpha=0.2, color=prop_cycle[0])
plot(gaussian_number_range, test_scores_mean[column],
'-', label = 'test',color=prop_cycle[1])
plt.fill_between(gaussian_number_range,
test_scores_mean[column] - test_scores_std[column],
test_scores_mean[column] + test_scores_std[column],
alpha=0.2,color=prop_cycle[1])
plt.xticks(gaussian_number_range)
print(column)
plt.locator_params(axis='y', nbins=5)
plt_configure(xlabel='Number of Normal Distribution', ylabel=column_name,
figsize=(3,2), legend={'loc':'best'})
plt.show()
In [ ]:
for fig in [fig_hist, fig_kde, fig_em, fig_gmm]:
display(fig)
for fig in [fig_time_variability_3d, fig_time_variability_cdf, fig_time_variability_weibull,
fig_adjecent_variability_3d, fig_adjecent_variability_cdf, fig_adjecent_variability_weibull,]:
display(fig)
In [ ]:
import time
save_notebook()
time.sleep(3)
location_name = get_location_name(ipath)
current_file = 'GMM.ipynb'
output_file = './output_HTML/'+location_name+'.html'
output_HTML(current_file, output_file)